clc;clear all;close all;
global alpha betta depr sigmaC phi govexp taumax N T ALPHA popweight
global kstst lstst cstst rstst wstst kstart kistart tauK0 mustst
global Delta1 Delta2 lamda transfer kg lam kmax
global restrict

%mainbenchNoL;

load('benchN2NoL.mat')

%find long-run steady state given initial guesses for lam, Delta1 and Delta2

kistart = [k1ss k2ss];
kstart = sum(popweight.*kistart)+kg;
kmax=5;

%ALPHA=1;

lamda=lam;

Delta1 = 0;
Delta2 = 0;

restrict=0
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

T=200; %we assume that the economy reaches the steady state after T periods

tauK0 = taumax;


%initial guesses 
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
    X(T+i) = mustst; %mu
    X(2*T+i) = max(1-0.1*i,0); %gamma
end;
transfer=4;
X(3*T+1:3*T+N) = [transfer lamda];
X0=X;

ALPHall=[2.9 3 3.08 3.09 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.55 3.6 3.62 3.625 3.63 3.631 3.632 3.65 3.7]

for ii=1:length(ALPHall)

ALPHA=ALPHall(ii)

if ii>1
    X=XX_fb(ii-1,:);
end
transfer = X(3*T+1);
lamda = X(3*T+2:3*T+N);
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
    
Delta1=0;
Delta2=Delta1;
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
X = fsolve('resid3vNoL',X,options)
  
X1=X;
X=X1;

if ii==1
    for i=1:T;
        time1 = i/T;
        X(i) = time1*(kstst) + (1-time1)*kstart; %k
    end
    X = fsolve('resid3vNoL',X,options)
end

XX_fb(ii,:)=X;

k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
transfer = X(3*T+1);
lamda = X(3*T+2:3*T+N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

tauh= 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*(kistart(1)*(1+(r(1)-depr)*(1-tauK0))-transfer))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];
plot(tauKt)

[swi dur]=min((tauKt-taumax/2).^2);
dur

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

if sum(V>=VsqNoL)==2
    PI(ii)=1
else
    PI(ii)=0
end
Vall_fb(ii,:)=V
VsqNoL
save('solN2NoL_fb.mat')

end

Vcap_fb=Vall_fb(:,2)
Vwor_fb=Vall_fb(:,1)
 
term1=(1-sigmaC)*(1-betta)*Vcap_fb;
epscap_fb=((term1+1).^(1/(1-sigmaC))/cistst_bm(2)-1)*100
term2=(1-sigmaC)*(1-betta)*Vwor_fb;
epswor_fb=((term2+1).^(1/(1-sigmaC))/cistst_bm(1)-1)*100

save('solN2NoL_fb.mat')


